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Abstract 

We propose a Langevin model with Coulomb friction. Through the analysis of the 
corresponding Fokker-Planck equation, we have obtained the steady velocity distri- 
bution function under the influence of the external field. 
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1 Introduction 



The Lan gevin equ ation is widely used to describe the motion of stochastic 
particles. ( Kubo92h In a typical situation a Brownian particle is driven by the 



random force originated from collisions of molecules in the environment and 
the friction force proportional to the velocity of the particle. 

The friction force is not always proportional to the velocity of particles. 
A well-known counter example can be found in the dry friction of macro- 
scopic materials, so called Coulomb friction in which the dynamical friction 
is independ ent of the speed of a material but depends on the direction of 
the motion. ( Perrson98h In general, the motion of such the material obeying 



Coulomb friction is not affected by the thermal fluctuations, because the mate- 
rial is too large to be fluctuated by the thermal force. However, if the material 
we consider is enough small like a nanomaterial or the thermal fluctuation is 
replaced by the mechanical fluctuation force, Coulomb friction force with fluc- 
tuating rando m force may b e relevant. In a recent paper, in fact, Kawarada 



and Havakawa (|Kawarada04l ) suggest that Langevin equation associated with 
Coulomb friction can be an effective equation in a vibrating system of granu- 
lar particles through their simulation of granular particles. They demonstrate 
that (i) granular particles in a quasi-two-dimensional container obey an ex- 
ponential velocity distribution function (VDF), (ii) the exponential VDF is 
produced by Coulomb slip between particles and fixed scatters because their 
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simulation without Coulomb slip produces a Gaussian-like VDF, (iii) Langevin 
equation with Coulomb friction produces an exponential VDF, and (iv) the 
relation between the diffusion coefficient D x and the granular temperature T 
in their simulation obeys D x oc VT as is expected from Langevin equation 
with Coulomb friction. Their result may encourage us to look for relevancy 
of such Langevin equation in physical systems, and investigate mathematical 
properties of the equation. 

In this paper, we investigate mathematical properties of Langevin equation 
with Coulomb friction. The organization of this paper is as follows. In the next 
section, we introduce our model. In section 3, we obtain the steady VDF under 
the influence of the steady external field. In section 4 we discuss our result 
through the comparison of our result with the simulation of driven granular 
particles in a quasi-two-dimensional box. In the final section, we conclude our 
results. 



2 Langevin model 



Let us consider the motion of a particle whose velocity is v subjected to the 
dynamical Coulomb friction force — £mv/|v| with the mass m. When we as- 
sume that the motion is affected by the external force F ex and the uncorrelated 
fast agitation by the environment, the equation of motion is given by 

<iv . v F ex . 
at |v] m 

where £ is the random force whose a component is assumed to satisfy 

< Ut) >= 0, < UtW) >= 2D8 a f8(t - t'). (2) 



Here < • • • > represents the average over the distribution of the random 
variable £ and the diffusion constant D in the velocity space satisfies the 
fluctuation-dissipation relation D = (y/T/(3m), where T is the temperature 
defined from the second moment of VDF as T(t) = \m J dvv 2 P(v, t). Since we 
assume that higher order correlations of £ vanish, Langevin equation (1) can 
be converted into Fokker-Planck equation for the probability function P(v, t): 

dP(v,t) F ex d n , . d ( v „. , „ d n , A 

+ • 7^( v > t) = -. (-P(v, t) + D-P(v, t) . (3) 



dt m dv dv \ Ivl dv 
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Here P(v, t) satisfies the normalization condition 

JdvP(v,t) = l. (4) 



Kawarada and Havakawa (|Kawarada04l ) indicate that eq.(3) without F ex has 



the steady solution obeying an exponential VDF. Namely, the VDF approaches 

rj 2 

P(v, t) -> <p (v) = — exp[-rjv] (5) 

Z7T 



as time goes on, where rj = (/D. Thus, our system has completely different 
properties from those in the conventional Langevin model. 

It is convenient to use dimensionless distribution function. For this purpose, 
we may introduce 

P(v, t) = v (ty d P(c, t), c = v/v (t) (6) 



with the normalization 

dcc 2 P (c,t) = 1 (7) 



where d is the dimension, and Pq is the scaled VDF without F ex . We adopt Pq 
instead of P in eq.(7), because the determination of the scaled factor is difficult 
if we use P. In fact, T is the function of F ex and cannot be determined without 
the complete form of VDF. 



3 Steady state distribution function 



Let us consider the situation that a particle is moving under the influence of 
F ex = mgz with the unit vector z parallel to the direction of the external force 
and Coulomb friction force associated with the random force. We assume that 
the system is a two-dimensional one, because the particle is typically located 
on a substrate when Coulomb friction is important. We are interested in the 
statistical steady state in the balance between the friction and the external 
force. 

Let 9 be the angle between the direction we consider and z, the steady equation 
becomes 



} dP + sin 2 6 dP _ P + dP + d 2 + 1 d + 1 d 2 

dv v d cos 9 v dv dv 2 v dv v 2 d6 2 
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where v — |v|. 

Now we assume the expansion 

oo 

P(y,9) = P n( v ) cosn9 = P o( v ) + p i(v) cos 9. (9) 

n=0 

Here the terms with n > 2 are irrelevant, because the contribution of n > 2 is 
orthogonal to those of n — 0, 1, and the effect of gravity appears in the term 
of n = 1. Thus, the normalization (4) is reduced to 

oo 

2tt J dvvP (v) = 1, (10) 





and eq.(8) becomes 



g [ cos 61* + cos 2 6P[ + —A] = C[— + P<^ + (— + A') cos 0] 

V V V 

+D[P" + ^ + (P" + S.) cose - ^ cos0], (11) 

t> t> V 

where P' n and P" represent dP n /dv and d 2 P n /dv 2 , respectively. From the in- 
tegrations of eq.(ll) multiplied by cosn9 with n — 0, 1 over (0, 27r), we obtain 



+ ^) = C( V + P o) + *W + f ), (12) 

^o = C(^ + A') + W + ^-§). (13) 

Equation (12) can be integrated easily as 

P = ^P + -Pi (14) 



for g 7^ 0. We can check that P oc e ^ and Pi = are the solution of (12) 
and (13) which correspond to the steady solution (5) at g — 0. 

To solve eq.(13) we adopt the form 

P (v) = f(v)e^ v (15) 
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and substitute this into (14), eq.(14) is reduced to 

Pi(v) = jf(v)e-*». (16) 

With the aid of eqs.(15) and (16), eq.(13) can be rewritten as 

/"' + ( I ~ v) f" " + f + = 0, (17) 

where e = (g/() 2 - 

Since we do not know the general procedure to obtain the solution of the 
third order ordinary differential equation (17) and it may not be easy to get 
the numerical solution for eq.(17) around the singular point v — 0, we rewrite 
it as a set of the second order differential equations under the assumption of 
small e. For this purpose, we expand 

f(v) = ^[l + ef^ + e'f^ + •••]. (18) 

Here the prefactor is determined by the normalization (10). Thus, we obtain 
the series of equations: 

h n 1 + (-- V )h' 1 -\h 1 = -ri 3 (19) 

v v 

K + (--v)h' n -^h n = -r ] 3 f^+r ] 2 h^ 1 for n > 2 (20) 

v v 

with h n (v) = df^ n \v)/dv. We also assume 

oo 

J dvvf {n) (v)e~ vv = for n > 1, (21) 
o 

which ensures that the normalization condition (10) is satisfied for any e. 

Equations (19) and (20) have a common homogeneous part supplemented by 
the inhomogeneous parts in the right hand sides. Thus, at first, let us obtain 
the solution of the homogeneous equation: 

w" + (-- nV - = 0. (22) 
v v 2 

The solution of eq.(22) can be obtained by the method of a series-expansion 
around the singular point v — 0. The result is 
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oo n 9 

^^(^r"^- 1 -" 1 (23) 



jj2 oo k -y 

w 2 (v) = —wi(v) \nv-J2 yiYl ~} v 

- k=l ' r=l 



n 2 e vv 

= — wAv) lnt> {7 + ln(r]v) — Ei(—r]v)} (24) 

2 v 

where 7 is Euler's constant 7 = lim n ^ oo(X^Li - — Inn) = 0.57721 • • • and 

00 . 
e 1 



Ei(-x) = - / dt— (25) 



is the exponent ial integra l function. Here we use Frobenius' method and 
Bessel's formula (jArfken95h 



00 k lx k 



E E ~TT = ^{7 + Inx - J5i(-x)}. (26) 

fe=l r=l 



The general solutions of inhomogeneous equation (19) or (20) is given by 
the sum of a special solution of the inhomogeneous equation and the linear 
combination of the homogeneous solutions (23) and (24). For 0(e) it is easy 
to confirm that the special solution h^s is given by 

hM = rfv. (27) 



On the other hand, the leading order contribution e vv /v should be canceled in 
the linear combination of W\ and w 2 to obtain the finite result for / °° dvvf( 1 >(v)e~' nv . 
Thus, we obtain the solution as 



hiiv) = cq{w 2 {v) + —{7 + h\rj)w 1 {v)} + rfv 

= c { Eli ~ Vv) - i±^( 7 + \n( V v))} + rfv (28) 
v v 

where c is a constant. We should note that hi(v) approaches hi(v) — > — CqT] 
in the limit of v — > 0. 

From the integration we obtain the first order solution as 

f {1) (v) = Co{GH (t]v\ q J Q ^ ~r]v(j-l) + (-f + r]v)\nr]v-^(\nr]v) 2 } 
+\v 2 + c i ( 29 ) 
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where c\ is a constant and 



dsz s n r(6 j - s) n r ( x - a 3 + s ) 

x f f[ r(i - b 3 + s) f[ r(o,- - s)] (30) 

Vj'=m+1 j=n+l / 



is Meijer's G function ([Meii er4 it ) which satisfies 



d 







-Ei(-z). 



(31) 



From the normalization condition (21) constants Co and c\ in eq.(29) become 

2 2 

c = 1, and ci = 4 - ^ + ^ - y ~ -2.98961, (32) 



respectively. Here we use 8\ = G^l ( 1 
formulae ( Wille86l ) 



-1 1 



23 v -i o 



1.64549.. from the 



f dmrp-^G 31 fnv\ 01 \ - r7 32 ( 1 

J dvv e Cr a3 I Q Q Q I - I/ C^ll, Q Q Q 







-a 1 



(33) 



Using this distribution function, the mobility \x defined through 

OO 2lT 



v = y dtw y d9v cos9P(v,9) = fiF e , 







(34) 



is given by 



T 



(35) 



where we use D = C y /T/3m (|Kawarada04l) . 



Now, we comment on the dimensionless forms of our result obtained here. 
With the aid of (6) and (7), it is easy to check -qv becomes y/Ec in the scaling 
form. In addition , the prefactor D/g in (16) becomes D/(gvo) = 1/y/Ee. 
Figure 1 plots the scaled Pq(c) and Pi(c) as functions of c = v/v . 
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P n (c) 




Fig. 1. The scaled Po(c) (solid line) and -Pi(c) (dashed line) as functions of c = v/vq, 
where P n (c) denotes -Po(c) or Pi(c) in the figure. For P\{c) we plot (c)e~^ c to 
remove the effect of the expansion parameter e. 

4 Discussion 



Let us discuss our result. First, we can obtain higher order terms such as f^(v) 
and f^(v), because the homogeneous solutions in eqs.(23) and (24) can be 
used in any order. However, these solutions may have complicated forms and 
we had better use numerical integrations to represent inhomogeneous terms. 
Because of the limitation of the length of this paper, we have omitted to give 
the explicit forms of higher order terms here. 



Second, Kawarada and Havakawa (|Kawarada04l ) assume that the diffusion co- 
efficient in the real space is proportional to the diffusion constant in the ve- 
locity space, but this assumption may not be true for our case. Let us briefly 
evaluate the diffusion constant in the real space as follows. The basic equation 
is 



d t P(r,Y,t) = TP=(T + T 1 )P 



(36) 



where 



T = DdJri- + d v ) and T ± = -v • V. 

v 



(37) 
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Now, we introduce the projection operator 

Ph(r, v, t) = (p (v) f dvh(r, v, t). (3£ 



With the aid of Mori-Zwanzig projection method) Zwanzig61; Mori65l ) . in gen- 
eral, the distribution function obeys an identity 



d t PP = PTPP + PT J dre {t - T)QT QTPP{T) 

o 

+PTe^ v QP{t ) (39) 

where Q = 1 — P. It is easy to show that the first and the last terms of right 
hand side of (39) are zero if the initial distribution is proportional to <po(v). 
We also assume that I\ can be regarded as the pertubation of r because 
the distribution function is unchanged within the mean-free path of particles. 
Thus, we may replace e^~ T ^ r by e^~ T ^ ro in the second term of right hand side 
of (39). Noting 

gr> (^(r, t) = -v • V<p {v)F(r } t) (40) 

v 

r v^o(^) = -C-v 2 o(«) (41) 

V 

with F = PP, eq.(39) may be reduced to 
t 

d t F = - J dvT J c/re ( ^ r)r °v^oVF 

o 

t 

~V J ' d\ J dre-v ( - t - T) v 2 (p VF 
o 

~ i J dvv 3 <p (v) V 2 F = D X V 2 F (42) 

in the long time limit t ^> i>q/C- Thus, the diffusion coefficient in the real space 
may be D x oc D 3 /^ oc T 3//2 . We should note that this procedure is not exact 
because t t> /C does not ensure the approximation /q dre~^ v ~ v/( for 
large v ^> v . However, the treatment may be plausible and is contradicted 
with the previous result of simulation of granular particles. 

We also note that the behavior of VDF under the influence of g in granu- 
lar particles does not coincide with the theoretical prediction presented here 
quantitatively. As shown in Fig.l, P and Pi have peaks around c = 1 and 
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become negative near c = 0. The behavior of the simulation for granular par- 
ticles is qualitatively similar, but the peak position is located near c = 0.5 and 
the negativity in the vicinity of c = is not large in the simulation. 

Thus, we will have to check whether the connection between Langevin equation 
with Coulomb friction and the vibrated granular particles confined in a quasi- 
two-dimensional box is superficial. The detailed quantitative comparisons will 
be reported elsewhere. 



5 Conclusion 

In conclusion, we have developed the theory of Langevin equation with Coulomb 
friction and obtained the steady solution of VDF under the influence of a 
steady external field. We also discuss the relation between our system with 
the vibrated granular systems confined in a quasi-two-dimensional box 

We would like to thank A. Kawarada for fruitful discussion. This work is par- 
tially supported by the Grant-in-Aid for Scientific Research (Grand No. 15540393) 
of the Ministry of Education, Culture, Sports, Science and Technology, Japan. 
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